Phase diagram of a two-dimensional lattice gas model of a ramp system 
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Using Monte Carlo Simulation and fundamental measure theory we study the phase diagram of a 
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I. INTRODUCTION 

The presence of liquid-liquid (LL) equilibrium in sim- 
ple fluids has drawn considerable attention in recent 
yearSji^ mainly due to its connection with the exis- 
tence of certain thermodynamic, structural and dynamic 
anomalies in liquid water The fact that there are sig- 
nificant regions of the phase diagram in which an increase 
of temperature at constant pressure is associated with a 
corresponding increase in density, or in which diffusivity 
is enhanced when the system is compressed, might be at 
first sight somewhat counterintuitive, and has therefore 
motivated a remarkable research effort. Most of the sys- 
tems that exhibit this peculiar behavior are also known 
to present low density solid phases (with coordination 
numbers ranging from 2 to 5) less dense than their liquid 
counterparts. Melting upon compression is a common 
feature that has to be accounted for as well. 

In this regard, simple models constitute a fundamen- 
tal aid that can allow to identify those essential features 
key to the presence of the aforementioned anomalous be- 
havior. Recent research has focused on two main cate- 
gories of models: orientational and isotropic. The former 
class of models is constructed bearing in mind the ori- 
entational character of the hydrogen bond interaction or 
the strong directional character of the covalcnt bonding 
characteristic in systems with low density solid phases, 
such as silica^ germanium oxide^ or phosphorus X For 
this class of materials, a series of realistic potentials have 
been employed in order to characterize their anomalies 
via computer simulationji^^iii^^ Useful as these stud- 
ies might be, a better insight can be gained from sim- 
pler models which can be dealt with in some cases even 
analytically. Perhaps the precursor of the simple ori- 
entational models is the Bell-Lavis two-dimensional lat- 
tice model of water^ recently somewhat extended by 
Barbosa and Henriques^ In addition to these, the Mer- 
cedes Benz model of water—, the 3D lattice gas model of 
Roberts and Debenedettii^, and the two-dimensional as- 
sociating lattice gas model of Hcnriques and Barbos a 15 i 16 



must also be mentioned. 

The complexity of the above mentioned orientational 
models can be further reduced. A weighted orientational 
average from these types of interactions would lead in 
most cases to isotropic models with several interaction 
ranges. And, since the pioneering work of Hemmer and 
Stelljil it turns out that the presence of two compet- 
ing scales or interaction ranges has been found to lie 
at the heart of the existence of multiple phase transi- 
tions in otherwise "simple" fluids . The ramp potential 
model proposed by Hemmer and Stell regained atten- 
tion when Jaglai^ stressed the similarities between its 
behavior and the anomalous properties of liquid water. 
Since then, a good number of works have been devoted to 
the continuous ramp potentiaL 19 i 20 ' 21 i 22 i 23 i 24 Other sim- 
ple models with competing ranges of interaction, such 
as the hard-sphere square shoulder-square well poten- 
tial have also been shown to exhibit LL equilibria. 25 
But not only continuous models can furnish an illustra- 
tive qualitative picture of the phase behavior and various 
anomalies found in water and related systems. Isotropic 
lattice gas models have proven to be able to describe 
the qualitative features of these systems rather accu- 
rately. One dimensional ) 26 ' 27 i 28 two dimensional 2 ^ and 
three dimensional 2 ^ lattice gas models have been stud- 
ied, using either mean field approaches, transfer matrix 
methods and/or computer simulation. 

In this paper we will consider a two-dimensional lat- 
tice gas model closely connected with the one studied in 
Ref. [28| in three dimensions. The model is characterized 
by two competing interaction ranges (a nearest neigh- 
bor hard core exclusion and a finite repulsive interaction 
on the next to nearest sites). This model is strongly re- 
lated to the continuum shoulder model studied in Ref. [2o1 . 
when the attractive interactions are absent. Our study 
will focus on the reentrant melting of the low density solid 
phase, using both computer simulation and Lattice Fun- 
damental Measure Theory (LFMT).^!^^^ We will 
see how the theoretical approach provides a qualitatively 
fairly approximate picture of our model phase diagram. 
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The rest of the paper is sketched as follows. A brief 
description of the model is introduced in the next section. 
Section III is devoted to the simulation methodology. De- 
tails on the finite size scaling analysis of the transitions 
are included in Section IV. Section V describes the LFMT 
as applied to this model, and finally our most significant 
results and conclusions are presented in Section VI. 

II. MODEL 

We consider a two-dimensional model defined on a tri- 
angular lattice. A given site of the lattice can be cither 
empty or occupied by one particle. The occupation of 
that site excludes the occupation of its six nearest neigh- 
bor (NN) sites. In addition there is a repulsive interaction 
between pairs of particles located in pairs of sites that are 
next to nearest neighbors (NNN). 

The potential energy of an acceptable configuration is 
then written as: 

U = e niTij. (1) 

<ij> 

where < ij > indicates the set of NNN pairs of sites; 
the coordinates arc equal to zero for empty sites and 
equal to one for occupied sites; and e > 0. 

In the limit of high temperature the interactions be- 
tween NNN become negligible and the system behaves as 
a hard core lattice gas with NN exclusion. Such a model 
is the well known hard hexagon model, which exhibits a 
continuous order-disorder transition. The location of this 
transition was obtained by Baxter t 35 i 36 and it is believed 
to belong to the same universality class as the 3-state 
Potts model in two dimensions^ At high density the 
system adopts an ordered structure (that will henceforth 
be referred to as T3) in which the sites occupy prefer- 
entially one of the three sublattices of the system [see 
Fig. (Ha)], with a triangular structure. The density at 
close packing is 1/3 (one third of the sites are occupied). 

In the low temperature and low pressure limit (equiv- 
alently e — > oo) another lattice hard core model is met. 
In this case an occupied site excludes the occupation of 
its NN and NNN. Under these exclusion rules, at high 
density we find again an ordered phase in which particles 
sit preferentially on one of the four corresponding sublat- 
tices [see Fig. QJb)]. The close packing density is in this 
case 1/4 (one fourth of the sites are occupied), and this 
ordered phase will henceforth be denoted T4. Attending 
to the symmetry of the order parameter and the dimen- 
sionality of the system, one expects to find an order- 
disorder transition belonging to the universality class of 
the 4-state Potts model in two dimensions.— Theoretical 
analysis 3 ^ and simulation result o 39 ' 40 have shown that 
this transition is also continuous. 

As in previous work by some of the authors^ we are 
interested in the phase transitions of the system at inter- 
mediate temperatures, where three phases: disordered, 
T4, and T3, can appear. 




FIG. 1: (a) Ordered structure at high temperature (large filled 
circles) and close packing. The three sublattices are identified 
by circles with different shading. Lines link nearest neighbor 
sites. The distance between nearest neighbors in each sublat- 
tice is \/3, i.e. the distance to next nearest neighbors in the 
original lattice, (b) Low density ordered structure appearing 
at low temperature (large filled circles). The four sublattices 
are identified by circles with different shading. Lines link 
nearest neighbor sites. The distance between nearest neigh- 
bors in each sublattices is 2, i.e. the distance to third neigh- 
bors in the original lattice. One fourth of the sites (large 
circles) of the lattice are occupied. 



III. SIMULATION METHODOLOGY 

In order to obtain the phase diagram of the system 
we have made use a number of Monte Carlo (MC) sim- 
ulation techniques. At high temperature we have used a 
flat- histogram algorithm ; 24 ' 2 ^ 41 ' 42 inspired on the Wang- 
Landau (WL) method , 43 ' 44 to compute the Helmholtz 
energy function for all possible densities of the systems 
at fixed temperature and volume. In practice, the simula- 
tion procedure can be regarded as an extension of the MC 
simulation in the grand-canonical ensemble (GCE), in 
which the different number of particles are not weighted 
by a fixed chemical potential, but using a weighting func- 
tion which permits us to extend the densities entering the 
sampling to an arbitrary range. Such a weighting func- 
tion is not known a priori, but it can be computed fol- 
lowing the strategies underlying the WL method. Once 
the Helmholtz energy function is known, for a given tem- 
perature, for all the densities and several system sizes, 
one can resort to finite-size scaling strategies to locate 
the phase transitions. This technique, with some modi- 
fications, has also been used to determine the transition 
occurring at low temperature and moderate density be- 
tween a fluid (F) disordered phase and a triangular T4 
phase. 
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In addition, we have made use of the so-called Gibbs- 
Duhem integration (GDI) techniqu e 45 ' 46 to determine the 
transition between the two ordered phases and to check 
the consistency of the results. 

In what follows we summarize these techniques. 



A. Flat-histogram simulation at constant 
temperature 

The flat histogram algorithm is divided in two 



parts r 



The first one is devoted to find a weighting 



function to sample efficiently a prefixed range of densi- 
ties, at constant temperature T, and volume. In this 
part we make use of the WL strategy. In the second part 
the actual sampling of the system properties is carried 
out. We will comment later about the specific details of 
each of these two steps. In the application of the algo- 
rithm, two types of moves, denoted as translation and 
insertion/deletion attempts, are considered. 

Translational moves are carried out as follows: (i) A 
particle is selected at random, and removed from its po- 
sition, Ri, in the system, (ii) A trial position R* (which 
could even be the previous one) is selected at random 
with equal probability from those positions which are nei- 
ther occupied nor excluded by the NN interaction, (iii) 
The new position is accepted with a probability given 
by the standard Metropolis criterion applied to the in- 
teractions of the particle in the current and in the trial 
positions, i.e 



45,47 



-4(R-|Ri) = min{ 1 



exp [-/%4(Rj)] 
exp[-/%(Ri)] 



(2) 



where lAi is the interaction energy of particle i, and (3 = 
1/fceT, with ks being Boltzmann's constant. 

In the second type of MC move, a change of the number 
of system particles, N, is attempted. First of all it is 
randomly decided (with equal probabilities) whether to 
increase or decrease N. Let us first consider the most 
common case in which N > 1 for the removal attempts 
and N < N max — 1 for the insertion attempts, A max being 
the maximum number of particles to be considered. If the 
number of particles is to be reduced, an occupied position 
is selected at random and its particle is either removed or 
left according to the acceptance criteria. If an insertion 
is attempted, a non-excluded position (if there is any, 
otherwise the insertion attempt is directly rejected) is 
selected to insert a particle, and as above the acceptance 
criteria are applied. 

In order to present the acceptance probabilities of these 
attempts in the context of a flat histogram procedure, 
let us first write the canonical configurational partition 
function, 



Q(N,M, T) = i ]T exp [-(3U(R N )] 



{R«} 



M 



N 



(3) 



JV! 



(exp [-/3U}) . 



where {R*} is the full set of M N possible configura- 
tions of TV distinguishable particles over a lattice with 
M positions. The factor M N JN\ is the contribution of 
the ideal lattice gas (system without interactions) and 
(exp [—(3U]) accounts for the excess contribution to the 
configuration integral. 

The sampling of different values of JV can be carried out 
by introducing a weight function uj(N). The probability 
of a given configuration of the system then becomes 



P(R N \M,T) cx w(JV)cxp [-PU{R N )] . 



(4) 



Integrating ^ over all the configurations of indistin- 
guishable particles for a given value of JV we get 



P(N\M, T) oc lu(N)Q{N, M, T). 



(5) 



For the particular choice u>{N) = z N , we obtain the prob- 
ability of JV in the GCE: 



P(N\fM, M, T) oc z N Q(N, M, T) 



(6) 



where z is the activity, which is related with the chemical 
potential, fx, as z = exp(/3fx). 

In order to perform an effective sampling of the ther- 
modynamics of a system for a wide range of densities, at 
fixed conditions of M and T, we can choose a weighting 
function different to that defining the GC ensemble. In 
practice we look for a prescription ujf(N) that produces 
a flat distribution P(N\M,T), i.e.: 



Lu f (N)(x 1/Q(N,M,T). 



(7) 



For practical purposes we introduce the function F(N) 



rN 



i.e. 



defined by w f {N) = JV! exp [.F(JV)] /M 
F ex (N,M,T)/k B T + K, with K being a constant, and 
F ex (N,M,T) the excess contribution to the Helmholtz 
energy function. 

We can now define the acceptance probabilities of an 
attempted change in the number of particles. For that 
we start off from the detailed balance equatio n 45 ' 47 ' 48 



P(R JV )VF(R iV+1 |R iV )^(R JV+i |R JV ) 

= P(R Ar+1 )VF(R Ar |R JV+1 )^(R JV |R Ar+1 ), 



(8) 



where and R w+1 are two configurations of the sys- 
tem with N and N+l particles, respectively, which differ 
just in the fact that the latter has an additional parti- 
cle with respect to the former. The right-hand side of 
([8]) is the probability of moving from the configuration 
R N to the configuration R*" 1 " 1 in a given Monte Carlo 
step, whereas the left-hand side corresponds to the prob- 
ability of the reverse move. By W(a|6) we denote the 
probability of choosing a as trial configuration in an in- 
sertion/deletion move, when the current configuration is 
b. Finally ^4(a|6) is the probability of accepting the MC 
move b — > a. 
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Taking into account the above described method of 
performing the insertion/deletion moves we can write 



W(K N+1 \R N ) 



W(R N \K N+1 ) 



1 



2A/- pos (R^r 
1 

2(7V+ 1)' 



(9) 
(10) 



where Af poa (R ) is the number of available positions 
(those that are not excluded due to hard core interac- 
tions) in the system with configuration R w . Taking into 
account ((TJ)— (jTDJ) we can write: 

■A(RA|R"+i) =«pH^+i)— ST" 

x exp[F(N + 1) -F(N)] 

where AUn+i is the change of the potential energy of the 
system when introducing the new particle. 

We have performed two types of calculation. At high 
temperature the full range of possible number of parti- 
cles, < N < M/3, has been sampled. In this case we 
have introduced transitions between the empty and the 
fully occupied lattice. This is feasible because the total 
number of configurations of both cases is known exactly, 
namely one for the empty lattice and three (correspond- 
ing to the filling of each of the tree sublattices) for the 
fully occupied lattice. Thus W{R°\R M/3 ) = 1/2 and 
VF(R°|R M / 3 ) = 1/6 and the corresponding acceptance 
ratios can be obtained for this particular case. We will 
call this scheme cyclic sampling. 

At low temperatures we found difficult to sample 
the whole range of densities because the WL procedure 
showed slow convergence. So in order to analyze the 
transition between the gas and the T4 phases we per- 
formed the simulations in the range < N < N/A. In 
this case a cyclic scheme is not feasible because at T > 
other configurations different to those of the perfect T4 
structure are possible for N = M/A. Therefore in the 
insertion/deletion sampling one directly rejects selected 
trials of particle deletion when N = and of particle 
insertion when N = M/A. 

Technical details about how to compute the Helmholtz 
energy fuction F(N,M,T) using flat-histogram tech- 
niques can be found elsewhere i 41 i 42 ' 49 Here we will just 
mention the basic ideas underlying the calculation. Sim- 
ulations are divided in two parts: equilibration and sam- 
pling. In the equilibration part F(N) is modified during 
the simulation run to push the system to visit all the 
values of N in the selected range. This equilibration is 
split in stages, each one run until certain convergence 
criteria are satisfied. As the stages go on the changes in 
^F(N) are smaller. At the end of the equilibration one ex- 
pects to have an appropriate estimation of !F(N). Once 
the equilibration part is finished, the resulting function 
J-(N) is kept fixed and the sampling part of the simula- 
tion starts. During this part one computes the probabil- 
ity of each value of N (from which a refined result for the 



Helmholtz energy function F(N, M, T) can be obtained) 
and different properties of the system, such as energy, 
energy fluctuation, order parameter, etc. The sampling 
part is divided into blocks in order to estimate error bars. 



B. Gibbs-Duhem integration 



As in previous work: 



,21. 2S 



with systems exhibiting a 



similar phase behavior, we have employed GD I 5 ' ° in 
the computation of the phase diagram. In the present 
case GDI was used to determine the T3-T4 transition 
and to check the consistency of the flat histogram MC 
calculations. 

The analysis of the Helmholtz energy function in the 
limit T — > leads to a value of p/e = 12 for the T3- 
T4 transition. After a number of short calculations we 
conclude that such a value hardly varies for low temper- 
atures. 

GDI goes as follows. For fixed volume (M) systems, 
the changes in the grand potential can be written as 



d(-(3 P M) = Ud{3 - Nd((3p). 



(12) 



At given (3 and p phase equilibrium exists if the pressure, 
p, is equal in both phases. Now if one changes, for in- 
stance, (3, the change in (3p to keep phase equilibrium is 
given by: 



. Ait ,„ 
Ap 



(13) 



where u = U/M, p = N/M is the density, and AX de- 
notes the the difference of the values of the property X 
in the two phases. Equation (| 13[) . or some variants of it, 
can be used to build up numerical integration schemes 
to compute the phase equilibrium of discontinuous tran- 
sitions. In practice we have performed GDI using two 
different integration schemes. At low temperature we 
have used 



dp 



Au 
~A P 



(14) 



whereas for the intermediate temperatures, like those at 
which the T3-T4 line is expected to meet the F-T4 line 
transforming it into a T3-F line, we employed 



dp = ^d((3p). 
Au 



(15) 



IV. ANALYSIS OF THE PHASE TRANSITIONS 

In order to locate the order-disorder phase transitions 
at high and low temperatures we can use the results for 
F(N, M, T) to obtain the probabilities in the GC ensem- 
ble, 

P{N\p, M, T) cx cxp {-/3F(N, M, T) + (3pN} . (16) 
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Then we search for the value of the chemical poten- 
tial, fj, c (M, T), that maximizes the density fluctuations. 
For this conditions we compute the average density, 
p c {M,T), and the momenta of the distribution of den- 
sities m n {M,T) = {(5p) n ) u ,T, with Sp = p- (p c (M,T)), 
for n = 2, 3, and 4. 

According to the definition of p c {M, T), we must have 
7713 (M, T) = 0. The system size dependence of rri2(M, T) 
and the ratio g 4 (A/,T) = m 4 (M, T)/[m 2 (M, T)] 2 , allow 
us to characterize the (possible) phase transition. 

In principle, given the symmetry of the model, one 
expects that at high temperature the phase transition 
will be continuous and belong to the universality class 
of the 3-state Potts model in two dimensions,— whereas 
at low temperature the order-disorder (F-T4) transition 
is expected to lie in the universality class of the 4-state 
Potts model in two dimensions. 

The scaling behavior that standard finite size scaling 
(FSS) predictsi£^£i goes as follows: 

Pc{L,T) ~ Alc (T)+a AlJ L- 1 ^, (17) 
p c {L,T) ~ Pc (T) + a p L 1 ^- d , (18) 
m 2 (L,T) ~ a m2 L a '»- d , (19) 

where we have used L (related with M by M = 2L 2 ) as 
the system length. These scaling laws are expected to be 
satisfied for large values of L. In these equations d is the 
dimension of the lattice (d = 2), and v and a are critical 
exponents, which are expected to take the values 3 ^ v = 
5/6, a = 1/3 for the F-T3 continuous transition, and 
v = 2/3, a = 2/3 for the F-T4 continuous transition. 

At intermediate temperatures the nature of the transi- 
tions can change, and eventually become first order. This 
fact can be studied by analyzing the behavior of g 4 (L, T) 
with the system size. In general, for discontinuous tran- 
sitions the value of (74, goes to 1 as L — > 00, signaling the 
presence of two well defined narrow peaks in the distri- 
bution probability of the density. 

V. THEORETICAL APPROACH 

We have performed a theoretical analysis of this model 
using LFMT.^ii 3 ^^ This theory is the lattice coun- 
terpart of Roscnfcld's Fundamental Measure Theory^ 2 - 
(FMT), and its construction is based on the approach 
through zero-dimensional (Od) cavities and dimensional 
crossovers of Tarazona and Rosenfeld.— In short this 
theory amounts to computing the exact functional for 
a certain set of small graphs (the Od cavities) and then 
build the simplest functional for the whole lattice which 
provides the exact result for density profiles that are 
zero everywhere in the lattice except in a Od cavity. In 
essence this theory is the grand-canonical functional ver- 
sion of Kikuchi's cluster variation method^ in Morita's 
formulation.— 

For short-range interacting lattice gases, the con- 
struction of a LFMT density functional is particularly 



simple^ 3 - Here we will explain in detail how to apply it 
to the concrete model we are studying. For a full account 
of the theory in all its details and with all its properties 
the reader is referred to Refs. [H andl34l. 

The starting point of the theory is the choice of a set of 
so-called maximal Od cavities. Zero-dimensional cavities 
are subgraphs of the lattice such that every two particles 
placed on them necessarily interact. They are "maximal" 
if adding a new node to the graph breaks down this Od 
requirement. If the interaction is purely hard-core exclu- 
sion, every Od cavity can hold just a single particle. If on 
top of that there is a soft interaction, then more than one 
particle can be present in a Od cavity. For the particular 
model we are considering, in a triangular lattice C, the 
set of maximal cavities is given by 

W 4 = |J W 4 (r), W 4 (r) = lj# 4). 

(20) 

The label r in the graphs denotes the position of the node 
beside it (the remaining nodes of the graph are labelled 
accordingly). Notice that cavities placed at different po- 
sitions in the triangular lattice are considered different. 

The second step is to complete this set by closing it 
with respect to non-empty intersections, i.e. if two over- 
lapping cavities are in the set, so must be their intersec- 
tion. The full set of cavities resulting from this operation, 
W, can be described as 

4 

W=|J>V l , (21) 

i=l 

where W4 is given by (f2"0"| , and similarly are defined VVj , 
i = 1, 2, 3, with 

W 3 (r) = { r &, y}, (22) 
W 2 (r) = { r c«,, y, (23) 

Wi(r) = {^}. F ' (24) 

The density-functional form that LFMT prescribes for 
this lattice gas is then F[p] = Fid[p] + F ex [p], where 

^idM = ^Kr)[mp(r)-l] (25) 

andH^ 

PF ex [p]= [-Mw(C,£)]$o(C), (26) 
cew 

with /zyy(C, C) a combinatorial object known as the 
Mobius function of the set W U {C}^^- and <J> (C) the 
excess free-energy density functional of the system when 
constrained to be within the cavity C. &o(C) is therefore 
a function of the density profile p(r) at the nodes of C 
alone. We will return to its calculation in brief. 

The Mobius coefficients yUw(C, C) satisfy the recursion 

ti w {C,L) = -I- Mw(C',£), (27) 

cccew 
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with which it can be obtained for any cavity C of the set 
W. It turns out that every cavity C £ is contained 



in the same number of cavities C £ W, with i < 



3- 



We shall denote this number Mij . Then pyy (C, £) is the 
same for all cavities of the same set Wi- So by denoting 



recursion (|27p becomes 

TO; = — 1 — Afjj TOj . 



(28) 



(29) 



It is easy to find that matrix M = (M^-) is 



M = 



/0 6 6 12\ 

2 5 

3 

\0 0, 



(30) 



On the other hand, it follows from (|2T[) that 7774 = — 1. 
This determines the remaining coefficients as m 3 = 2. 
m-2 = 0, and mi = — 1, and therefore the functional as 

[3F C M= E $ o(C)-2 E *o(C) + E $ "( C )- 

cew 4 cew 3 ceWi 

(31) 

Let us now compute the functions $o(C) for all C £ W. 
Actually it is enough to obtain this function only for 
the maximal cavities (those of W4), because any other 
cavity C must be — by construction — a subgraph of one 
of the maximal cavities, and therefore $o(C) can just be 
obtained by setting p(r) = at the nodes of the maximal 
cavity which do not belong to C. On the other hand, 
by symmetry the functional dependence of <&o (C) on the 
densities at the nodes of C will be the same for all the 
maximal cavities of this model. In other words, the only 
function we need to obtain is 



(32) 



where the cavity nodes are labelled generically and 
should be appropriately replaced by the nodes of the 
corresponding cavity. Hence this is a function of p = 

(pi, P2, Pa, Pi)- 

We start off by writing down the grand-canonical par- 
tition function for such a cavity, namely 



5 = 1 + z + KZ1Z4, 



(33) 



where k = 
any e > 0), 



o /3[M-Kxt(0] 



e ^ e (therefore we have < k < 1 for 
z = Zi denotes the total activity, and 
Vext(i) representing any external field 
acting on node i. In obtaining (f3"3")l we have made use of 
the fact that the cavity can accommodate at most two 
particles, and this only if they occupy nodes 1 and 4. In 
that case they interact through the soft potential. From 
(|33p we can obtain the densities as pi = (zi/ T E,)d'E,/dzi, 



and the correlation between nodes 1 and 4 as P14 
(ziz i /E)d 2 E/dzidz4. This yields 



^1(4) — z l(4) + KZiZ 4 , 
S P2(3) = z 2(3), 
Spi4 = KZ1Z4. 



(34) 
(35) 
(36) 



Adding the equations for pi up we obtain 

Sp = z + 2kziZ4 = 5 — 1 + Spi4, (37) 
where p = YliPi- Then 

- = l-p + p M . (38) 

On the other hand, from Eqs. (f3~4"|) (|36|) it follows that 
£1(4) = S(pi(4) — P14), hence substituting this expressions 
in ([36]) and using |38|) we obtain 

k(Pi - pu){pi - pu) = Pia{1 - P + Pu)- (39) 
The solution to this second-order equation for P14 is 

Pu(p) = 2( 1 ^ K ) { ~ 1 + P- K (Pi + Pi) 

+V / P ~P + «(Pi + Pi)} 2 + 4k(1 - K) Pl pij , 

(40) 

Finally we get $0 through a Legendre transform, i.e. 



*o(p) = E P* ln ( z */Pi) ~ ln3 



i=i 



:p+(l-p) ln(l -P + P14) 

i=l,4 V Pl J 



(41) 



As explained above, we obtain the $0 for the non- 
maximal cavities by setting p,; = at the corresponding 
nodes. This leads to 

$ (C) =p+(l-p)ln(l-p), p = £>, (42) 

for all CeWiUW 2 U W 3 . 

Equations (25|), ^ and gDJ-gl]) complete the pre- 
scription for the density functional. 



VI. RESULTS 

A. Monte Carlo Simulation results 

The simulations were carried out on rectangular lat- 
tices of different sizes, built by replicating in both direc- 
tions those shown in Fig. Q] which depict rectangular 
lattices containing M = 2 x L x L sites (with L = 6). 



7 



0.00 
0.10 
0.20 
0.30 
0.40 
0.50 
0.63 



2.406(2) 
3.152(2) 
3.927(2) 
4.728(3) 
5.546(9) 
6.393(3) 
7.510(11) 



Vc 

0.8279(9) 

0.8430(10) 

0.8561(11) 

0.8671(6) 

0.8765(20) 

0.8840(6) 

0.892(4) 



TABLE I: Computed points for the F-T3 transition. Error 
bars are given in brackets, in units of the last figure of the 
property, and correspond to a confidence level of about 95%. 



1. F-T3 transition 

For the location of the F-T3 transition we used the WL 
cyclic sampling described in Sec. MI A| for values of fie = 
0.00, 0.10, 0.20, 0.30, 0.40, 0.50, and 0.63. For each of 
these values simulations were carried out for ten system 
sizes: L = 12, 18, 24, 30, 36, 42, 48, 54, and 60. We com- 
puted the pseudo-critical quantities fi/j, c (L,T), p c (L,T), 
etc, and extrapolated these data to the thermodynamic 
limit. To take into account possible deviations from the 
scaling laws due to the relatively small system sizes, we 
use the ad-hoc fitting 



X C (L,T) =X C (T) 



fc=i 



(43) 



where X c represent some physical property at the tran- 
sition point, and b x the critical exponent appearing in 
its corresponding scaling law [c.f. Eqs. (fT7l) - p8|) ]. The 
number m is chosen to be either one or two, according to 
a chi-square test<££ 

We have found that for most values of fie the pseudo- 
critical chemical potentials /j, c (L,T) can be fitted for the 
whole set of system sizes using a second-degree poly- 
nomial [m = 2 in Eq. ([43]) ]. In the particular case of 
fie = 0.63 the smallest system sizes (L = 12, L = 18) 
were discarded due to the interference with the F-T4 
transitions. 

The results for the F-T3 transition are collected in 
Table [J (notice that the densities are expressed in terms 
of packing fractions, i.e. 77 = 3(N)/M), and plotted in 
Figs. 2(a)| [3j and [4] The particular case fie = corre- 
sponds to the so-called hard hexagon model, whose crit- 
ical properties are known exactlyi 35 i 36 The exact values 
are f3fi c = log [(11 + 5\/5)/2] w 2.4061, and ? ]c = 3p c = 
3(5 — v5)/10 ~ 0.8292. A comparison with the extrapo- 
lations given in Table |T] suggests that the latter are quite 
accurate — though not perfect — in the estimation of the 
exact values. 

The results for r) c (/3e) and fip c (fie) are well represented 
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FIG. 2: Temperature-density phase diagram from Monte 
Carlo simulations (a) and LFMT (b). The three phases are 
labeled F for the fluid, and T3 and T4 for the two solids. 
Coexistence regions are marked with two labels. 



by the fits 

?7F-T3(/3e) 
/3^ F -T3(/3e) 



0.8280 + 0.1582/3e - 0.0922(/3e) 2 ,(44) 
2.4056 + 7.3074/3e + 1.6079(/3e) 2 
-0.5478(/3e) 3 . (45) 

2. F-T4 transition 



As in the previous case, we have selected a number of 
representative temperatures and performed simulations 
for several system sizes. In most cases we considered the 
same sizes (L=12, 18, ■ • • , 60) as for the F-T3 transition. 
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FIG. 3: Pressure-temperature phase diagram. The vertical 
axis represents the pressure divided by fcsT (ftp) and the hor- 
izontal axis represents the temperature in the form e _/3l! , with 
e > the soft repulsion between NNN sites. Symbols are the 
Monte Carlo simulations and lines represent the predictions 
of the LFMT. The shaded region is a set of values that cannot 
be reached because the pressure jumps discontinuously at the 
T3-T4 transition. 




FIG. 4: Same as Fig. [3] but for the chemical potential fx. 
The shaded region is a set of values that cannot be reached 
because the chemical potential jumps discontinuously at the 
T3-T4 transition. 



f3e 




Vc 


1000 


1.756(1) 


0.560(2) 


10 


1.756(2) 


0.560(2) 


5 


1.774(1) 


0.562(2) 


4 


1.806(1) 


0.565(2) 


3 


1.897(1) 


0.573(2) 


2 


2.165(1) 


0.596(2) 


1.5 


2.491(2) 


0.619(2) 


1.25 


2.777(3) 


0.637(4) 


1.0 


3.246(4) 


0.658(4) 


0.9 


3.541(3) 


0.676(5) 


0.8 


3.950(2) 


0.688(4) 


0.75 


4.230(2) 


0.695(3) 



TABLE II: Computed points for the continuous F-T4 tran- 
sition. Error bars are given in brackets, in units of the last 
figure of the property, and correspond to a confidence level of 
about 95%. 



In addition to /3/z G (L,T) and p c (L,T) we also payed at- 
tention to the quantities m 2 (L,T) and g±{L,T) because 
this transition appears to be discontinuous in some cases. 

The precise location of the multi- critical point 
— where the transition changes from continuous to 
discontinuous — is a hard task due to the first order tran- 
sition being quite weak. In order to make an approximate 
estimation we focused on the changes of gi{L, T) with L 
and found that the change of character of the transition 
occurs at about f3e « 0.75 ± 0.05. Surprisingly, it is pre- 
cisely at this temperature that the simulation results are 
well represented by the scaling law (fT9|) . 

For /3e > 0.75 we have estimated the critical proper- 
ties of the F-T4 line in the thermodynamic limit using 
the same strategy as for the F-T3 case, employing the 
critical exponents of the 4-state Potts model in two di- 
mensions. The results are gathered in Table [TTJ The 
results at low temperature show a good agreement with 
those reported by Zhang and Dcng^S: /3/i c = 1.75682(2), 
and Tjc = 0.540(12). 

For /3e < 0.75 the plot of the chemical potential as 
a function of rj shows a loop for the different system 
sizes considered. This is a signature of a first order 
phase transition. The change in density is quite small, 
and the transition is rather weak. Taking this into ac- 
count we have estimated the location of the transitions 
by fitting the results for each simulated temperature 
to equations of the form with b x = 1, m = 2. 

The properties considered were fi c (L,T), r] c (L,T), and 
Ar] c (L,T) = y/ m>z(L, T). Notice that, in general, for 
discontinuous transitions lim£-»oo m,2(L, T) ^ 0. Thus, 
in the thermodynamic limit, the packing fractions of the 
two coexisting phases are then obtained as 



Vf(T) = Vc (T) - A Vc (T) 
VTi{T) = T] C (T) + A?] C (T) 



(46) 
(47) 



The results for the discontinuous F-T4 transition are 
collected in Table HIU 
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pe 


0.70 


0.65 


0.63 


P/1F-T4 


4.61(2) 


5.18(2) 


5.59(4) 


Vf 


0.705(2) 


0.717(2) 


0.724(2) 


VT4 


0.718(2) 


0.733(2) 


0.737(2) 



TABLE III: Computed points for the discontinuous F-T4 
transition. Error bars are given between parentheses, in units 
of the last figure of the property and correspond to a confi- 
dence level of about 95 % 



B. Gibbs-Duhem Integration 

After a number of tests we analyzed the discontinuous 
T3-T4 transitions and its continuation as F-T4 transi- 
tion using GDI with lattices of size L = 120. This rela- 
tively large size was chosen since the end of this line at 
low values of (3p corresponds to the equilibrium between 
a low density fluid and the T4 phase. This transition 
is relatively weak and because of this we found that for 
smaller values of L one of the subsystems often under- 
goes a phase transition and both subsystems end up in 
the same phase, with the corresponding breakdown of 
the GDI scheme. The use of large simulation boxes then 
makes possible us to reach values of (ip that allowed us 
to check the consistency between GDI and WL simula- 
tions. Technical details of the integration algorithm can 
be found elsewhere i 24 ' 28 The integration steps and length 
of the simulations were chosen after performing a number 
of tests. 

The GDI was divided in two parts. In the first part 
we perform an integration following the scheme given in 
Eq. |H1). The first point was chosen to be T* = (/3e)- 1 = 
0.35, fi/e = 12. The integration was carried out using 
a temperature step of AT* = 0.025. At each step we 
run long simulations (2 x 10 6 cycles, each cycle including 
M/3 insertion/deletion attempts; averages are taken over 
the second half of the simulation). The integration was 
carried out up to T* = 1.25 (or /3e = 0.80), where we get 
coexistence for p/e = 11.947 ± 0.001. 

The second part of the integration was carried out us- 
ing the scheme given in Eq. (fT5|) . The starting point 
was that defined by the previous integration (T* = 1.25, 
j3e = 0.80, ftp = 9.558). The integration step was then 
A((3p) = —0.020, and the length of the simulations was 
about 1 x 10 6 cycles (for each system and step, averag- 
ing over the second half of each run). As in the previous 
part, a number of additional simulations with larger in- 
tegration steps, smaller system sizes, and shorter runs, 
were carried out in order to test the integration accuracy 
and to estimate error bars. 

Simulations were launched to execute 201 integration 
steps (to reach, in principle, a final value of f3p = 5.558). 
We observed that for this line GDI required both large 
systems and precise estimates of the integrand in order 
to avoid the collapse of the method before reaching the 
F-T4 discontinuous equilibrium found with WL Simula- 



Point 


TP (T4-F-T3) 


HT (F-T4) 


DC (F-T4) 


pe 


0.660(2) 


0.619(1) 


~ 0.75 


Pn 


7.78(2) 


6.35(2) 


~ 4.0 




0.892(5) 


0.745(1) 


~ 0.69 




0.892(5) 






T)TA 


0.7491(1) 


0.745(1) 


~ 0.69 



TABLE IV: Singular points of the phase diagram: TP stands 
for triple point, HT for high temperature end-point, and DC 
for the change from discontinuous to continuous behavior of 
the fluid- T4 transition. Error bars are estimated by compar- 
ing results from different GDI trajectories. 



tions. For instance, using large simulations (precise inte- 
grands) with L = 60, A.(j3fJt) = —0.05, it was possible to 
get good estimates both for the triple point T4-F-T3 and 
for the temperature at which the density of phases F and 
T3 are equal at equilibrium, but shortly after reaching 
the latter point the algorithm failed due to the transi- 
tion of one of the phases into the other. With L = 120 
the integration stayed stable until reaching the expected 
final value of ftp, ~ 5.56, which was found to be con- 
sistent, within statistical uncertainty, with the value of 
the F-T4 equilibrium obtained through WL simulation 
at f3e = 0.63. 

From the results obtained with GDI, together with 
those previously reported on the F-T3 transition, we can 
estimate the position of two of the special points in the 
phase diagram, namely the triple point T4-F-T3 and the 
point of maximum temperature for the F-T4 equilibrium 
(at this point both phases have the same density). These 
results, together with those of the change of the transition 
order of the F-T4 equilibrium, are collected in Table [TV] 



C. LFMT calculations 

We shall now apply the density functional obtained 
in Sec. [V] to determine the temperature-density phase 
diagram of this fluid. To this purpose we have to study 
the three phases involved: the uniform fluid and the two 
solid lattices (T3 and T4, c.f. Fig. [J). 



1. Uniform fluid 

If every site has the same average occupancy p (hence a 
packing fraction rj = 3p), the free energy per unit volume 
(in fcgT units) will be 

$ =p\np+ (1 - p) ln(l -p) - 4(1 - 3p) ln(l - 3p) 

+ 3(1 - 4p) ln(l - 4p + p u ) + 6p In (\ - ^ , 

(48) 
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where p u = p u (p, P, P, p) is 

P 14 = 2^(2(2-^-1 

- 4(2 - k) P + 4(4 - 3k)p 2 } 
The pressure can then be obtained as 



(49) 



PP = P' 



dp 



= In 



{\-ipY 



(l-Ap + p u ) 3 (l-p) 



(50) 



and the chemical potential is given by (3p = ($ + /3p)/p. 



2. T3 solid phase 

The T3 solid occupies one of the three sublattice shown 
in Fig. QJa). F rom this figure it is clear that the den- 
sity profile will take one out of two values: The sites 
of the occupied sublattice — marked with large circles in 
Fig.rija) — will have a higher density pa, whereas the re- 
maining ones will have a lower density ps (by symmetry 
this value is the same for all those sites). The packing 
fraction will then be 77 = pa + 2ps- 

If we regard the lattice as a tiling of rhombii, we can 
realize that there arc three different kinds of them: Those 



with pa at positions 1 and 4 and ps at positions 2 and 3, 
those with pa only at position 2, and those with pa only 
at position 3. There is the same number of each kind. 
The contribution of the two latter to the rhombic cavities 
will be the same, but different from the contribution of 
the former. Thus 



go(g) 
V 



X] ~ TT ' = ®o{pa, Pb, Pb, pa)+2<P {pb, Pa, Pb, Pb), 
cew 4 

(51) 

with $o(p) given by (|4Tj) . As for the triangular cavities, 
J2i Pi — Pa + %Pb = 77 for any of them, so 

-^^ = 2?/ + 2(1 -77) 111(1-77). (52) 



cew 3 



V 



Dimcrs do not contribute to the functional, so the last 
contribution will be 



2^ x t = ~ + oC 1 ~ P^)ln(l - pa) 



V 



3 3 



+ 3(1- Pb) Ml- Pb)- 



(53) 



Putting all together and adding the ideal part the re- 
sult is 



Hpb; v) =-(i-v + 2p B ) Mi -v + 2 pb) + -(1 - pb) Mi - pb) - 4(1 - v) Mi - v) 

+ (1 - 2/y + 2p B ) ln(l - 277 + 2p B + Ci) + 2(1 - 77 - p B ) ln(l - v -p B + &) (54) 
+ (77 - 2 PB ) 1 2 111(77 - 2pB - 6) - § HV ~ 2Pb)| + 2p B Uln(p B - &) - | Inp B } , 

I 



where we have eliminated pa — 77 — 2p B , and 

Ci = Pu(pa,Pb,Pb,Pa), (55) 
6 = pu(pb,pa,pb,pb)- (56) 

The free energy as a function of 77 is obtained by mini- 
mizing the function (|54|) with respect to ps (always tak- 
ing the solution ps < ??/3). For k > n t = 0.5086 . . . (see 
below) there is a first order transition from a homoge- 
neous fluid to the T3 solid [see Fig. |2(b)| . In the limit 
K = 1 where the model becomes identical to the hard 
hexagon model, we find a wide first order transition. This 
model has bee n p reviously solved within the LFMT ap- 
proach in Ref. [32. LFMT is a mean-field-like theory, and 
therefore all second order transitions have a parabolic 
behavior of the order parameter — corresponding to a 
critical exponent (3 = 1/2. The exponent of the hard 
hexagon model (like that of the 3-state Potts model) 
i o 35 ' 36 (3 = 1/9. Looking at Figure 7 of Ref. HH one must 
admit that a discontinuous function is a better approx- 



imation to the behavior of the order parameter than a 
parabolic, mean-field one. So although not quite satisfy- 
ing, the result is quantitatively not too inaccurate. This 
also reflects in the fact that the pressure and chemical 
potential at the transition is rather close to the exact 
value, and it remains so for all n t < k < 1, as Figs. [3] and 
|4] show. 



3. T4 solid phase 

The T4 solid occupies one of the four sublattice shown 
in Fig.QJb). Again the density profile will take either the 
value pa at the sites of the occupied sublattice — marked 
with large circles in Fig. QJb) — or the value ps at the 
remaining sites (the same for all of them, by symmetry). 
The packing fraction will now be 77 = pa + 3p B . 

As a tiling of rhombii the lattice contains four kinds of 
them, each with an A site at one of the four positions. 
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FIG. 5: Free energy (in ksT units) per unit volume as a func- 
tion of the packing fraction 77, for re = e^ 13 '' = 0.3679, below 
the triple point. Dotted line represents the free energy of the 
uniform fluid; solid line that of the T4 phase, and dashed line 
that of the T3 phase. The filled square represents the bifur- 
cation point F-T3, and the filled circle the bifurcation point 
T3-T4. Notice the discontinuity of the derivative of the free 
energy at 77 = 0.75, the close packing of the NNN exclusion 
lattice gas. For r\ > 0.75 the free energy of the T4 phase is 
concave, so T4-T3 coexistence always occurs with a T4 phase 
at 77 = 0.75, and both the chemical potential and the pressure 
jump discontinuously at this transition. 



There is the same amount of each type. Thus 



$o(C) 3 
2^ — 77— = 7;®o(Pa,Pb,Pb,Pb) 



cew 4 



V 



(57) 



2®o(pb,Pa,Pb,Pb)- 



As for the triangles, one fourth of them have an A site 
and two B sites, and three fourths have three B sites, so 



\- $o(C) PA + 2PB , l n 



2pi 



cew 3 



x ln(l - pa- 2p B ) + -p B 
+ |(l-3p B )ln(l-3p B ). 



(58) 



Finally, the contribution of the point-like cavities 
(those of Wi) is 



PA) 



ceWi 



(59) 



+ -(l-p B )ln(l-p B ), 



because 1/4 of the sites are of type A and 3/4 of type B. 

Putting all together and adding the ideal part the re- 
sult is 



<£(p B ;r/) = ^-^ 



+ Ai + In 1 



4r? 
~3~ 



A- 



___ln( y -3p fl -A 1 



^-{Hp B - Ai) + 2 HPB - A 2 )} + Pb) ln(l - p B ) - (1 - 3p B ) ln(l - 3p B ) 



15 



p B ln p B - 3 1 



4// 

y 



PB In 1 



4// 



1 



1 



4// 



3pi 



ln 1 



477 

y 



1 (y - 3pB ) 1x1 (y - 3pB 



3pi 



(60) 



where we have eliminated pa — rj — 3ps, and 

Ai = pi4,(pa,pb,pb,pb), (61) 

A2 = Pl4 (pB , PA , PB , PB ) • (62) 

As for T3, the free energy of the equilibrium phase is 
obtained by minimization of this function with respect 
to pb (always choosing the solution ps < f?/4). In the 
limit k = the model is equivalent to a lattice gas with 
NNN exclusion. In this limit we obtain a wide first or- 
der transition from a uniform fluid to a T4 solid, which 
again is found to be continuous in the simulations. The 



values of the pressure and chemical potential for this 
transition are nevertheless rather accurately predicted 
(see Figs. [3] and HJ, so the same considerations as for 
the F-T3 transition in the hard hexagon (re = 1) limit 
hold here. We find a first-order F-T4 transition all the 
way up to k c = 0.5403 . . . , where it coalesces to an end- 
point (at 77 c = 0.7405 . . . ). It is to be noticed that this 
point is obtained with high accuracy (simulations yield 
k c fa 0.538 and 7/ c w 0.745; see Table HV]) . and that sim- 
ulations also find a first-order F-T4 transition near this 
point (see Fig. [J). 
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FIG. 6: Temperature-density phase diagram for the 3d coun- 
terpart of the model discussed in this article, on a simple cubic 
lattice. 

When k = the model has a close-packing at rj = 3/4; 
however, for any n > this limit can be crossed, although 
at a very high energetic cost. Once this cost is paid, the 
system greatly diminishes its entropy by reordering itself 
in a T3 structure. Hence the transition that is found, for 
all k up to not too far from n tl between a close-packed 
T4 solid and a nearly closed-packed T3 solid, which is 
also observed in the simulations (see Fig. [2]) ■ What hap- 
pens with the free energy at this T4-T3 transition is very 
peculiar and it is illustrated in Fig. [5l The derivative of 
the free energy per unit volume (in fcgT units) with re- 
spect to the packing fraction is discontinuous at 77 = 3/4. 
The free energy of the T4 phase is concave beyond this 
point, so there is a T4-T3 coexistence, but it does not 
satisfy the standard conditions of phase equilibria. As 
a matter of fact, both the pressure and the chemical po- 
tential jump discontinuously at this point. This creates a 
"forbidden" region in the pressure-temperature and the 
chemical potential-temperature, which can be observed 
in Figs. [3] and H 

At the value Kt there is a F-T3-T4 triple point, which 
the theory predicts very close to the value obtained in 
the simulations, n t sa 0.517 (see Table |TV|) . 

In the range n t < k < k c a reentrant T4-F transition is 
found before the F-T3 transition occurs. This reentrant 
behavior also appears in the simulations, although the 
coexistence region is wider because the F-T3 transition 
is continuous (see Fig. ■ 



VII. DISCUSSION 

From Figs. [2j El and 21 we see that the theoretical 
results for the equilibrium between the ordered phases 
are in reasonably good agreement with the simulation. 
The theoretical estimations of the triple point and the 
high temperature end-point of the F-T4 equilibrium are 
also well described. On the other hand the LFMT de- 



scription of the order-disorder transitions, both at low 
and high temperature is less accurate. The theory pre- 
dicts first order transitions whereas they are continuous 
in both limits. We have argued that this is so because 
the behavior of the order parameter is very sharp (the 
critical exponent f3 = 1/9), so much that a discontin- 
uous function is a better approximation to it than the 
simple parabolic behavior predicted by any mean-field- 
like theory (like this one). The theory can in principle 
be refined by considering larger "maximal cavities" . Al- 
though it would produce a far more complicated theory 
than the one presented here, and it is expected that its 
accuracy would increase, it is doubtful that the order of 
the transition would be corrected. No matter how much 
we complicate the theory, it does not cure its mean-field 
behavior near the transitions. This is also the reason why 
first-order transitions are described much better. 

In favor of this argument is the fact that, overall, the 
agreement between simulation and theoretical transition 
lines in the planes P/J.—T and @p-T is rather good. 

When comparing the phase diagram for the two- 
dimensional system on the triangular lattice with that 
of the three-dimensional systemsS [c.f. Fig. [6] in a simple 
cubic lattice we find a couple of qualitative differences. 
The first one concerns the nature of the order-disorder 
transition at low temperature, which is continuous in 
two dimensions and discontinuous in three dimensions. 
This difference can be understood in terms of the differ- 
ent dimensionality, and can be related with the critical 
behavior of Potts models^ in two and three dimensions. 

The second relevant difference arises when comparing 
the transitions between the two ordered phases. In two 
dimensions, at solid-solid equilibrium the high density 
solid is nearly close packed for a wide range of temper- 
atures, whereas this is not the case for its three dimen- 
sional counterpart even at the lowest temperatures. This 
difference can be explained as follows. At low temper- 
ature phase equilibrium is essentially controlled by the 
condition of minimum energy. The energy per unit vol- 
ume of the closed packed configuration in two dimen- 
sions is u* = U/Me = 3. At slighly lower 77 there are 
vacancies. The way in which they minimize the energy 
is by not being nearest neighbors in the solid lattice (i.e. 
next-nearest neighbors in the underlying lattice). This 
way each vacancy reduces the energy by 6e. Thus, for 
2/3 < r\ < 1, u*(rj) = 2r\ — 1. If we consider, at the same 
density, a system separated into a close-packed T3 and a 
T4 phase, then N T3 + A T4 = N and 3iV T3 + 47V T4 = M, 
and the energy of this system will be U = 3eiVx3. Hence 
u*(t)) = irj - 3. Since 4?y - 3 < 2rj - 1 for all 77 < 1, 
then the phase separated system is energetically favored. 
The same argument for the three-dimensional model of 
H0ye et al£2- yields the same energy, u* (77) ~ 677 — 3 (for 
3/4 < 77 < 1) in both cases, so in the three dimensional 
system the entropy does play a role in defining the den- 
sity of the high density solid, and the number of vacancies 
does not go to zero when approaching T = 0. 

The similarity between the phase diagrams in two and 
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three dimensions is remarkable, though. Peculiar fea- 
tures like the reentrant fluid phase or the vertical line 
at the closest packing of the loose solid when it coexists 
with the dense one, appear in both cases. 

Much to our surprise, we have realized that the LFMT 
for the three-dimensional model is far more complicated 
than that of the two-dimensional one. The reason is that 
the soft repulsion at NNN allows for maximal cavities 
with up to four particles at the same time. This not 
only introduces a much larger set of cavities to elaborate 
the density functional, but also the corresponding expres- 
sions for the $o functions is very cumbersome and hard 
to handle. On the other hand, given the similarity be- 
tween the phase diagrams, it seems that the physics of the 
model is already well captured by the two-dimensional 



version. 
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